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Abstract 

We obtained the spectrum of the Sinai billiard as the zeroes of a secular 
equation, which is based on the scattering matrix of a related scattering 
problem. We show that this quantization method provides an efficient 
numerical scheme, and its implementation for the present case gives a few 
thousands of levels without encountering any serious difficulty. We use 
the numerical data to check some approximations which are essential for 
the derivation of a semiclassical quantization method based also on this 
scattering approach. 


1 Introduction 

The Sinai Billiard (SB) has become a benchmark for testing ideas and methods 
in “quantum chaos”. Bohigas and his coworkers [1] calculated the lowest few 
hundreds of levels by using the boundary integral method [2], and showed that 
random matrix theory accounts very well for the numerical level statistics. Berry 
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[3] used the KKR method [4] as a starting point for his investigation of the 
semiclassical limit. More recently, the scattering approach to quantization [5] 
was demonstrated by applying it for the quantization of the SB. Last, but not 
least, this system was the hrst to be investigated experimentally by measuring // 
- wave reflection from a billiard shaped cavity [7]. 

In the present paper we would like to present yet another study of the SB, 
aiming this time to check and demonstrate the scattering approach to quantiza¬ 
tion in greater detail and depth. The scattering approach to quantization is based 
on the observation [5] that every billiard interior problem can be viewed upon as 
a scattering (exterior) problem, and the spectrum can be uniquely deduced from 
the knowledge of the corresponding scattering operator. As a matter of fact, one 
can show that the spectrum of the billiard coincides with the set of energies for 
which an eigenvalue of the scattering operator takes the value 1. 

The scattering system associated with a given billiard may be constructed 
in various ways. One may consider the billiard boundary as an obstacle in a 
scattering (exterior) problem [5]. Dirichlet (or Neumann) boundary conditions 
dehne uniquely the scattering operator and through Billet’s theorem [6], the in¬ 
terior energies are given as the energies at which an eigenvalue of the scattering 
operator obtains the value 1. The other approach, which will be pursued and 
investigated here, is based upon attaching appropriate “channels” to the billiard 
and the scattering operator coincides with the transmission or reflection matrix, 
corresponding to this quasi f-d scattering problem [5]. As a matter of fact, the 
exact quantization involves the extension of the concept of the scattering matrix 
to the space of “closed” or “evanescent” modes. When this is done one can prove 
(see chapter (2) that the spectrum of the interior problem coincides with the 
energies, for which an eigenvalue of the extended scattering operator takes the 
value 1. 

The scattering approach has a few advantages which will be demonstrated 
in the present note. One certain advantage is the numerical efficiency of this 
method. In chapter (3) we shall discuss some numerical results, and use the 
calculated spectrum and S matrices to demonstrate a few important relations. 
The scattering approach can be used naturally as a basis for the derivation of 
the semiclassical limit. To get to this limit, one has to go through two stages. In 
the hrst (the “semiquantal” approximation [5]) one neglects all the contributions 
from evanescent modes. The scattering operator is reduced to a hnite matrix and 
the secular equation for the interior problem becomes: 

Z,,{E) = det{I-S{E)) = {) (1) 

Here, S{E) is a A X A unitary matrix, whose dimension A is in general a step 
function of the energy, increasing by one at each threshold energy when a “closed” 
mode starts to conduct. The accuracy of the semiquantal approximation is ex¬ 
pected to deteriorate at the vicinity of threshold energies [5], in a way which will 
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be discussed in chapter (4). The semiclassical approximation is introduced by 
evaluating (1) semiclassically. One can show (see chapter (4)) that (1) can be 
expressed in terms of Tr[S^)^ 1 < < A, and the traces Tr[S^) can be written 

in terms of periodic orbits of the billiard! One derives in this fashion a semiclas¬ 
sical secular equation based on a hnite number of periodic orbits, without having 
to worry about divergent series and the proliferation of the number of periodic 
orbits. 

We shall end this paper with a summary, in which we shall list some inter¬ 
esting open problems and possible directions where the scattering approach to 
quantization might compete advantageously with other methods. We shall also 
discuss the relationships between the scattering approach and the semiclassical 
methods of Bogomolny [8] and others. 

2 Exact Quantization 

We consider a quarter of Sinai’s billiard, enclosed by the four sides of the square 
with diagonal {0,0} - {L^ L} and a quarter circle of radius R < L centered at 
the origin (cmp. Fig. 1). We are looking for a wave function vanishing 

at the boundary and satisfying 

(A + A;2)T^ = 0 (2) 



Figure 1: The Sinai billiard and the attached channel. 
In the region x > Rii can be decomposed into normal modes 
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The wave function will then take the form 
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where and are some constants. The total energy splits into two parts 






nvr 2 


L 


(5) 


where kn is the wave number for the motion in x - direction. Depending on n it 
will be a real or purely imaginary number. For n < A(A;) with 


A(A;) = [kLlir] 


( 6 ) 


we have = \kn\ and the mode is called open (traveling) whereas for a closed 
(evanescent) mode with n > A(A;) kn obeys k'^ = +z|A;„|. The number of 
open channels A(A;) depends on the total energy. The open modes consist of an 
incoming and an outgoing plane wave, the closed ones contain an exponentially 
increasing and a decaying contribution. For the sake of a compact formulation 
the terms incoming and outgoing will also be used for the evanescent modes. 

To be an eigenfunction of the closed billiard, the wave function must vanish 
ai X = L and this yields a system of linear equations for the coefficients bn'. 
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It will have a solution if 


detiSmn - e^^^-^Skn) = 0 . ( 8 ) 

In order to use this equation for the evaluation of the billiard eigenenergies one 
must know the S^n ^ function of k. They can be dehned as the coefficients 
of the outgoing modes m in the unique solution with a single normalized 
incoming mode in channel n 

^-iknX ^ikmX 

V m V ^rri 

In the limit x —)• oo and for n < A[k) the only terms that survive in (9) are 
the open modes. Hence S^n^ < A{k) is the unitary scattering matrix. The 

unrestricted S^n i® ^ generalization of the scattering operator for the scattering 
from the left billiard wall [21]. Obviously is a generalization 

of the scattering matrix for scattering from the right (straight) billiard wall at 
X = L in the same sense, and so we can write 

det{I - S^S^) = ^. (10) 

This is the exact secular equation for the billiard eigenvalues in the form of (1). 
For a usual scattering matrix unitarity and symmetry are important features 
and provide e. g. a valuable means for checking the numerical accuracy of the 
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computations. In appendix B we show how these relations can be extended to 
hold for the generalized S as well. 

It is well known [9] that the elements of the unitary S - matrix are related to 
the matrix elements of the Green function in a mixed coordinate - momentum 
representation. If we impose outgoing boundary conditions, the Green function 
G'“’“ takes the asymptotic form: 


G+{x,y;x\y) = - 
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ikmX+ik I x' \ 


( 11 ) 


where gmm’ are constants and depend on an additional boundary or symmetry 
condition imposed on G~^. For example, gmm’ = ^mm' will provide Gg , the Green 
function for the “pure” channel closed at x = 0 with a straight wall. In our case 
gmm’ are such that the Green function obeys Dirichlet boundary conditions on 
the left closure of the channel (x = 0, j/ > i?) and x^ -\- y^ = . Comparing (11) 

and (9) and using Green’s theorem we obtain: 


^mn g^ 


( 12 ) 


This can be also expressed as: 

(I-S’-U,, = X (13) 

I 

[ [ dydy'(f)m{y){G+{L,y,L,i/)-G^{L,y,L,i/))(f)m’{y') 

Jo Jo 

Another straight forward application of Green theorem yields: 

G+(f ",r') - G+(f ",r') = y G+(f",r)^G+(f,r')d 7 (f) (14) 

where 7 is the quarter of the circle = i?^, which is traversed in the 

positive mathematical sense. Equation (14) can be solved formally by successive 
iterations, resulting in: 


G+(f, ro) - G'^{f,fo) 


00 „ „ 

^ / d7(fi)-- - / d7(f)v)G+(f,f)v) 
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G^ (rjv_A', rjv_A'_i ) 


jv=i 

iV-l Q 


n 

K=0 


dvN-K 


(15) 


Equations (15) and (14) form a convenient starting point for the semiclassical 
approximation of . 

We shall now derive an explicit expression for the S matrix using the methods 
which were developed in [10]. There, the scattering matrix for an inhnite one¬ 
dimensional array of non overlapping spherically symmetric scatterers located at 
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n_De^ in the two dimensional plane is calculated (cmp. Fig. 2). Our generalized S' 
will be obtained by using the hard disk scattering phase shifts, and by extending 
the method to evanescent modes followed by a proper desymmetrization to fulhll 
the boundary conditions at the channel walls. Geometry requires for this purpose 
D = 2L. 


y 





Figure 2: 


The array of scatterers and the desymmetrized section of interest. 


Scattering off the inhnite periodic array of discs picks certain Bragg manifolds 
a with discrete y - momenta 


27r 

kl = Q + n —. 


(16) 


The boundary conditions will be fulhlled, if our wave function is antisymmetric 
with respect to x = 0 and both, J/ = 0 and y = i)/2. Within one Bragg manifold 
this can only be achieved for a = 0. The wave function is then periodic with 
the period D of the lattice and T(a;,0) = 0 will imply T(a;,i)/2) = 0. In the 
following we restrict ourselves to the case a = 0. 

A wave function with an incident plane wave from x = —oo in 

channel n is given in [10]-(II.14, 15) in normal mode representation which is valid 
for \x\ > R. In our notation it reads 


= iA “P (‘TT!/) + exp {,^y) (x < i?) 

>1':. = E {x>R} (17) 


m = — oo 


In [10]-(IV.5) the angular momentum decomposition 


y) = J2 aniRi{r)R‘'*’ (18) 

i 


is introduced which is valid for r < i)/2. The radial wave function is to vanish at 
the edge of the discs and is therefore determined by the T - matrix of the discs 
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which is diagonal in angular momentum representation: 

R,Xr) = J,{kr) + T,H+{kr]. 


( 19 ) 


Ji and H'^ are Bessel and Hankel functions, respectively [11], Writing the Lipp- 
mann - Schwinger equation for the problem in cylindrical coordinates it is possible 
to derive a system of linear equations determining all the coefficients We will 
not repeat this here but refer to [10]-(IV.9). However, in our case, where we want 
to include evanescent modes n is no more restricted by energy conservation but 
runs over all integer values. 

The only new ingredient required by this is a general angular momentum 
decomposition formula for exponentials: 

+ 00 

exp(zA;^a; + ikyy) = ^ (20) 

l= — oo 

Signal) ikx + ky 

k 

which is valid for both real and imaginary k^ and ky as long as k'^ = k^ -\- ky (see 
appendix A). 

For the coefficients in the decomposition of exp[—iknX + i^y) we will use the 
symbol Rni and from (21) we have 

R.nl = {-lyRn-l. ( 22 ) 

Proceeding in the same way which led to [10]-(IV.9) we now obtain 

dni’ {Sw — Fi_ii{kD)Tii} = Rn-i (23) 

i' 

with 

a-ni = {-lyan-i (24) 

apparent from symmetry. The structure function Fl[x) comes from the angular 

momentum decomposition of the Lattice Green function 

OO 

9k^\^'F') = X] G\y\f-,r'-\-nDFy) (0 < |r'| < |r| < D/2) (25) 

n= — oo 

where 

p) ^ iTrHy^\k\f-P\) (26) 

is the Green function of the plane. The structure function is dehned through 

OO 7 

Fi{x) = 2 ^ cos{—)H+{nx) (27) 

n=l 
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which is zero for odd 1. The evaluation of this structure function is the most 
difficult part of the applied method. It is adapted from [10] where it is described 
in appendix B. 

We take now antisymmetric combinations of the in order to match our 
boundary conditions and assume in the following x > 0: 

(28) 


^n{x, y) = (-x, y) - (x, y) (29) 

The normalization was chosen to yield (9) for the normal mode decomposition 
after using the symmetry of the problem to restrict the sum over modes to positive 
integers n, m = 1, 2, .... Now the S - matrix we are interested in is connected 
to the reflection and transmission amplitudes of the inhnite array by 


sL 




mn i — mn 


mn ^—mn 


(30) 


The angular representation turns into 


'^n{x,y) 



- an-i)Ri{r) sin Icf) 

i 


(31) 


where / = 2, 4, ... is now an index running over positive even angular momenta. 
From the linear system (23) we derive using the symmetries with respect to 
inverted indices 


— (In-l')iSu' + (Rl+l' ~ Rl-l')Tl') — Rn-l — Rnl ■ (32) 

I 


The angular momentum representation of the wave function in the vicinity of the 
discs can now be obtained by solving this linear set of equations. However, in 
order to get the S - matrix we need the normal mode representation of the wave 
function. For this purpose we use the fact that any two solutions of (2) obey 



for an arbitrary area Q. We choose Q to be enclosed by the channel wall y = L, 
the line x = L and the circle around the origin with radius L (cmp. Fig. 3). For 
the wave functions we take and T 2 = sinkm’X sin They both are 

antisymmetric with respect to j/ = T, therefore this line does not contribute to 
the integral and we obtain ly = 0 with 



( 34 ) 


8 




Figure 3: The desymmetrized section, the area Q is hatched. 



From these two formulas we extract 



(35) 

(36) 


For a compact notation we introduce the matrices: 



2i 

^nl — /j-— 

VknT> 

(37) 


2i 

Cnl — /7 - —{Rnl Rn—l) 

y'kniJ 

(38) 


dll' 

Piv = 7^ - 1“ Pi+i' ~ Pi-i’ 

J-v 

(39) 

(37) and (36) take now 

the forms 



AP = -C 

(40) 


= I-CA^ 

(41) 

Putting them together 

we end up with our hnal expression for 



= I + CP-^C^. 

(42) 


The above relations give explicit expressions for in terms of known quan¬ 
tities. They form the basis for the numerical evaluation of . 
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3 Numerical Results 


We start the presentation of our numerical results with an investigation of the 
accuracy of the eigenvalues we calculated. For this purpose we employ the sym¬ 
metry of the billiard with respect to reflections from the diagonal x = y. In the 
scattering approach this symmetry is broken by attaching the channel to one of 
the the billiard walls. However, the eigenvalues ki of the billiard must belong 
to one of the two representations of the symmetry. So they can be divided into 
two subgroups kf and k~ with the wave function inside the billiard symmetric 
or antisymmetric, respectively: 

{x^y) = {y-,^) (43) 

In polar coordinates the symmetric wave functions will contain only angular mo¬ 
menta I = i = 1,2,... If we denote the projection onto this subspace by Q'^ 
and the projection onto the complementary set of angular momenta by Q~ we 
have therefore in terms of the matrices (37) 

Aikt) = Q+Aikt) 0 = Q-Aikt) (44) 

Using S - matrices calculated with restricted space of angular momenta 

= I- CQ^A'^ (45) 

we have thus at the symmetric eigenvalues kf' besides 

Z{kt) := det(/ - S^S^ik+y) = 0 (46) 

also 

Z+{k+) := det (/-5'^5'^+(A:+)) = 0 (47) 

Though we cannot make use of the symmetry for a simplihcation of the calculation 
of eigenvalues (e. g. by searching only for all the k~ as in [1] and [3]), we obtain 
at least a means for estimating the deviations of the eigenvalues k^ computed 
by searching the zeros of (46) from the correct values k^: The numerical errors 
will cause deviations of Z“’“(A;4) from zero and the corresponding shift in the 
eigenvalue reads 

Sk=\kt-kt\ ~ (48) 

dk^ \kt 

Obviously, analogous reasoning holds for the k~. 

It is appropriate to measure the error estimate in units of the mean level 
spacing around the considered eigenvalue. The circles in hg. (4) show this relative 
error in a logarithmic plot for two different radii of the circle^. The error is about 
10“® for a quite small circle with R = 0.4 and about 10“® for a very large circle 


10 




Figure 4: Estimate of the numerical error at some eigenvalues for (a) R = 0.4 
and (b) R = 0.9 The eigenvalues were obtained including f closed channels. 
The circles show the estimated deviation in the wave number Ski normalized to 
the mean level spacing. The crosses denote the corre.sponding violation of the 
symmetry relations for the generalized S. 

(R = 0.9). We conclude that the scattering approach is well suited for an accurate 
calculation of thousands of eigenvalues for any given radius. The crosses in hg. 
(4) demonstrate the relative error in the generalized symmetry and unitarity 
relations for at the eigenvalues. They were obtained by picking the maximum 
relative deviation of a single matrix element of from the prediction of (80) - 
(82) and (89). We see that the violation of these relations is about 10“^^ and 
thus much smaller than the deviations in the eigenvalues. This means, that the 
most important error is introduced through the imaginary part of the structure 
function (27) which does not influence the checked relations (compare appendix 

B). 

The relative deviations of the eigenvalues are expected to increase with grow¬ 
ing radius of the circle and growing energy, since both will cause more angular 
momenta to contribute to the decomposition of the wave function around the 
circle: 

Imax ~ Rk (49) 

and the structure function in turn is needed up to 2lmax- Thus the maximum 
wave number k^ax where the relative deviation of eigenvalues is well below the 
mean level spacing decreases with growing radius of the disc. The number of 
eigenvalues up to k^ax decreases even faster, since growing radius of the inside 
circle means shrinking area of the billiard. For example, at i? = 0.4 it is possible 
to obtain about 5,000 eigenvalues for each of the two symmetry classes whereas 
for R = 0.75 we are restricted to the lowest 3,000. 

^In this and in all the following plots the edge length of the square was hxed at 4 = 1. 


If 





Once we have established the accuracy of the method, we use the resulting 
data to check various statistical properties of the spectrum. 



i 

Figure 5: Deviation of the eigenvalues from Weyl’s law. N(ki) — i at 200 eigen¬ 
values for R = 0.75. The dashed line is a running average over 200 eigenvalues. 


The spectrum of eigenvalues ki of the billiard can be described by the spectral 
density 

OO 

dfk) = '^6{k - ki). (50) 

8 = 1 

The integral of the spectral density 

k 

N{k) = J dk' dfk') (51) 

0 


increases by one at each eigenvalue and is called spectral staircase. The mean 
density of the eigenvalue spectrum in an interval containing many levels obeys a 
generalization of Weyl’s law [14] which reads in terms of the spectral staircase 


N{k) 



u 

dvr 


k + C. 


1,52) 


A is the area of the billiard, u the length of the circumference and C a constant 
which can be obtained from geometry as well. The oscillations of the actual 
spectral staircase N{k) around this mean are shown in hg. (5) for a fraction of 
the spectrum containing 200 eigenvalues. The mean value of 200 neighboring 
eigenvalues is plotted with a dashed line and is almost constant zero. 

The oscillations of N{k) — N[k) can be described semiclassically using the 
periodic orbits of the system. Besides unstable and isolated orbits the SB also 
permits orbits which are neutral and form continuous families. These so called 
bouncing ball orbits are not reflected from the quarter circle but only from the 
straight walls of the billiard. In [15] it is shown how the contributions of both 
unstable and bouncing ball orbits can be used to obtain a semiclassical prediction 
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for the length spectrum which is a hnite Fourier transform of the oscillating part 
of the spectral density: 


D{x, kr, 


dk coskx \d{k) — d{k) 


(53) 


This function is plotted in hg. (6) using 3350 eigenvalues for a SB with R = 0.75 
The lengths of the shortest periodic orbits are shown in hg. (6) by vertical 
lines. Each of these lengths corresponds to a peak in the spectrum. For the 
given radius, the only possible family of bouncing ball orbits is formed by orbits 
rehected back and forth between opposite sides of the billiard. These orbits have 
the length x = 2 and cause the most prominent peak in the length spectrum. A 
semiclassical prediction for the length spectrum was obtained using input from 
just the shortest ten periodic orbits and is displayed with a dashed line. It is a 
good approximation to the Fourier transform of the eigenvalue density. 



Figure 6: Length spectrum obtained from 3350 eigenvalues for R = 0.75 and 
kmax = 200 (full line) versus a semiclassical prediction containing the shortest 10 
periodic orbits of the billiard (dashed line). The vertical lines denote the lengths 
of the used periodic orbits. 


The SB is known to exhibit classically chaotic dynamics for any positive value 
of the inside radius R. Therefore the eigenvalue spectrum of the billiard can be 
described by the Gaussian orthogonal ensemble of random matrices GOE. To 
check this we employ the statistics of spacings between neighboring eigenvalues 
and the number variance S(iv) and rigidity A3(iv) which contain the two point 
correlations in the spectrum. The spectrum is unfolded using the generalized 
Weyl’s law and the statistics are compared to the predictions of random matrix 
theory([16], [17]). 
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Figure 7: 1600 antisymmetric eigenvalues for R = 0.75 (a) Nearest neighbor 
spacinq distribution (b) integrated distribution. The smooth line is the prediction 
from GOE. 

Fig. (7a, b) shows that the nearest neighbor spacings of the eigenvalue spec¬ 
trum are in good agreement with the GOF. The two point measures, however, 
follow the random matrix prediction only up to a maximum length and saturate 
above it. This behavior has been explained by Berry [18] using semiclassical 
arguments (see hg. (8)). 

In [10] and [12] it was shown that a S - matrix corresponding to a classically 
chaotic scattering system with time reversal symmetry can well be described by 
the circular orthogonal ensemble (COF) of random matrices. We would like to 
check some of the predictions of random matrix theory with the unitary part of 
the S - matrix employed for the billiard quantization. 

Consider hrst the distribution of the eigenphases of on the unit circle. For 
COF we expect a constant density On the other hand approaches / as the 
the radius shrinks to zero and for any radius we have a contribution to S from 
direct reactions. Therefore we have a peak in the eigenphase density around zero 
which falls in magnitude for increasing radius (hg. 9). We are not interested, 
however, in the contributions from direct reactions and will use therefore 
where they are not present. Fig. (10) shows that the corresponding eigenphase 
distribution is constant 

Ztt 

The distribution of spacings between neighboring eigenphases is expected to 
undergo a transition from an exponential behavior p[s) = to the Wigner 
surmise p[.s) = ^se~~^ as the dynamics of the corresponding classical system 
changes from regular to chaotic. Indeed, hg. (11) shows a distribution close to 
the Wigner surmise for large radius and a Poisson - like behavior for small discs. 

Two point correlations in the spectrum of eigenphases provide another test of 
the random matrix hypothesis for for the S - matrix. They can be expressed in 
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Figure 8: 1600 antisymmetric eigenvalues for R = 0.75 (a) Number variance (b) 
As statistics. The dashed - dotted line denotes the GOE prediction, the dashed 
line corresponds to a Poisson process. 
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Figure 9: Density distribution of the eigenphases of for (a) R = 0.4, (b) 
R = 0.65, (c) R = 0.9 Energy average over fOO matrices with 50 open channels. 


terms of the probability density P 2 {g) of finding two eigenphases with a distance 
g. The Fourier coefficients of this density are related to the eigenvalues 
[12, 13] through 

277 

5 ?^) := (A-l) J dgP2{g)E'^^ 

0 

= ^l^rh'^l^-l (54) 

A semiclassical description of the behavior of the coefficients in agreement 
with the COF prediction [16] has been derived in [12]. Apart from a non generic 
region for very small N one expects ~ ^ ~ 1 < A < A and ~ 0 

for N A. From hg. (12a - c) one can see that this prediction describes the S 
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Figure 10: Density distribution of the eigenphases of for (a) R = 0.4, (b) 

R = 0.65, (c) R = 0.9 Energy average over fOO matrices with 50 open channels. 
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Figure 11: Nearest neighbor spacing distribution of the eigenphases of for 

(a) R = 0.4, (b) R = 0.65, (c) R = 0.9 Energy average over fOO matrices with 
50 open channels. 


- matrix of the SB with large radius very well. For small radius the S - matrix 
can not be described by random matrix theory. 


4 Semiclassical Quantization 

As was explained in the introduction, the semiclassical approximation is derived 
in two steps. In the hrst, we restrict the S matrix in the secular equation (10) 
to the space of open channels. This is the “semiquantal” approximation for the 
secular equation, since the restricted S matrix is the fully quantal S matrix which 
describes the reflection as measured at large distances. This restriction is a natu¬ 
ral step within the semiclassical theory, because evanescent modes correspond to 
propagation with imaginary momenta, and generically their contribution is expo¬ 
nentially small in h. However, in the present case, the imaginary momenta start 
at zero at threshold, so that for any hnite value of h, there exists a corresponding 
neighborhood of the threshold energies, where the decay of the evanescent mode 
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Figure 12: Fourier coefficients of the cluster function as a function of N. (a) 
R = 0.4, (b) R = 0.65, (c) R = 0.9 Energy average over fOO matrices with 50 
open channels. The dashed line shows the prediction of the COE. 


which is about to be opened is too slow to justify its elimination. To check the 
accuracy of the semiquantal approximation, we used the methods and numerical 
data which were discussed in Chapter 3. above. 



47 50 53 


k 

Figure 13: E.stimation of the numerical error Ski for R = 0.75 The eigenvalues 
were obtained including 2 closed modes (full line) and no closed modes (squares). 
The crosses show the difference between both sets of eigenvalues. All data are in 
units of the mean level spacing. 

In hg. (13) we plotted all eigenvalues between the wave numbers where the 
15th and the 17th channel open. These thresholds are marked with vertical 
dashed lines. The full line shows the estimated deviation between the eigenvalues 
calculated with two closed channels and the exact locations of eigenvalues. It 
has been checked that inclusion of more evanescent modes does not alter the 
numerical eigenvalues for the parameters under consideration but neglecting these 
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two lowest closed modes does: The squares show the deviations of the eigenvalues 
obtained with unitary S from those with 2 closed channels and the crosses mark 
the corresponding error estimate according to (48). Both are of the same order 
of magnitude as it should be if (48) is correct. We see that the error coming from 
neglecting evanescent modes is in general not very large. It rises exponentially 
with decreasing distance to the next threshold where a new channel starts to 
conduct as predicted in [5]. 

The same behavior is also apparent from hg. (14). Here, the estimated 
relative error is shown for two different radii. The dashed line corresponds to 
the values obtained from 2 closed channels, squares and crosses mark the results 
for 0 and 1 closed channels, respectively. For the lower radius R = 0.75 there is 
no difference between the values for 1 and 2 closed channels but for large radius 
higher evanescent modes can also contribute as one can see from hg. (14b). For 
large radius the inhuence of the evanescent modes is in general bigger. This can 
be understood since the evanescent modes decay exponentially in the channel 
region R < x < L and should contribute less when the distance for the decay 
increases. In [5] the situation is different. There the channel was attached at the 
opening between the circle and the billiard wall. If the radius grows the distance 
for the decay increases and the inhuence of evanescent modes must diminish. 



Figure 14: Estimation of the numerical error Ski for (a) R = 0.75, (b) R = 
0.9 The eigenvalues were computed with 0 (squares), 1 (crosses) and 2 (line) 
evanescent modes and the error is in units of the mean level spacing. 


In [5] it is proposed that the error from neglecting closed channels should scale 
with the normalized absolute value of the wave number of the hrst evanescent 
mode 

A _ \kA+i\L 

7rV2ATT ^ ^ 


which is a number between 0 and 1. In hg. (15) we show the deviations of 
1000 eigenvalues computed without closed channels from those obtained with 
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Figure 15: Deviations between eigenvalues obtained with no and with one closed 
channel. First 1000 eigenvalues of the billiard with R = 0.75 The errors are 
normalized to the mean level spacing and the abscissa is the magnitude of k^+i 
in units of Tr/d{2A + 


one closed channel versus 6. The error increases exponentially as the magnitude 
of the lowest closed mode falls to zero. The error will be of the order of the 
mean level spacing if one eigenvalue is very close below a channel opening. The 
eigenvalues are, however, concentrated in regions where the error is very small. 

To end the discussion of the influence of evanescent modes we show the be¬ 
havior of the eigenvalues Si of in the vicinity of a channel opening. In hg. 

(16a) the magnitude of all the eigenvalues of is shown around the threshold 

for the opening of the second channel which is marked with a vertical line. Below 
the threshold there is only one open channel and therefore just one eigenvalue 
of is on the unit circle. Far away from the threshold the magnitude of all 

the other eigenvalues is very small and therefore they can be neglected. When k 
approaches the threshold the magnitude of the lowest evanescent mode increases 
from 0 to 1. For the same energy region in hg. (16b) the behavior of the eigen¬ 
values is shown in the complex plane. All the eigenvalues corresponding to closed 
channels are located in the origin. There is one eigenvalue corresponding to the 
hrst open channel traveling at constant speed on the unit circle. The eigenvalue 
corresponding to the second open channel tends to 1 as the energy approaches 
the threshold. This means, that the secular equation (1) will have additional so¬ 
lutions at the channel openings which do not correspond to billiard eigenenergies. 
The derivative of the eigenvalue of the opening channel is discontinuous at the 
threshold. 
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Figure 16: Behavior of the eigenvalues of for R = 0.75 in the vicinity of 

the opening of the second channel (vertical line), (a) magnitude (h) location in 
the complex plane 


Now that we have assessed the quality of the semiquantal approximation, we 
shall proceed to discuss some of its properties. The semiquantal secular equation 
reads 

Zsq{k) = det(/ — S{k)) = 0, (56) 

where the dimension of S[k) is A{k) = [^]. It can be written as 

Z., = jZiW (57) 

1=0 

in terms of fi{k)^ the coefficients of the characteristic polynomial of S{k). The 
unitarity of S{k) implies the important symmetry 

exp[-z0(A;)/2]/; = exp[z0(A;)/2]/;_; (58) 

where Q[k) = Ofk) — Att, and are the eigenvalues of S[k). The quan¬ 

tity Q{k)!{2 tt) can be shown to give the smooth number function of the spec¬ 
trum. This follows from the observation that for k values away from thresholds, 
the spectral density is given by 

A. 1 90(‘) , 1 A 1 a3rr(S"(n)) 

L + - L -- ^ - 

1=1 n=l 

The hrst term on the rhs is the smooth part of the density. Another important 
message of this expression is that the fluctuating part of the spectral density 
depends exclusively on Tr[S^{k)). The important role of TrS^ is even better 
emphasized when Newton’s identities are used to express the coefficients /;, 

TrA™ + fA-iTrS'^-^ + ... + fA-^+iTrS + mfA-m = 0; m = 1,..., A. (60) 
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We conclude, therefore, that the secular equation and hence the spectrum are 
determined once TrS^{k) are given for f < n < [(A + l)/2]. The second step in 
the derivation of the semiclassical quantization is to express these quantities in 
terms of classical periodic orbits. 

To get the semiclassical expression for TrS^{k)^ we shall hrst calculate Sl in 
the semiclassical approximation, the starting point being equation (14) of chapter 
(2). Consider hrst the free Green function 


G^{xy,x'y' 


^ <i>'m{y)<i>'m{y ) ^^^km\x-x'\ _ ^lkm(x + x')'^ 
2-1 k 

m = l 


( 61 ) 


Applying Poisson summation and performing the integrals by the stationary 
phase approximation, one gets the semiclassical expression 


^ (87r)i/2 ^ (AiC(f,r 


i[kls{r,r 


( 62 ) 


where s is an index which counts the classical paths from r to r ' which undergo 
on their way rehections from the boundaries at x = 0 , j/ = 0 and y = L. 
Note that the number of rehections from the line x = 0 can be at most one. 
Is is the length of the path. It should be noted that there exist also complex 
solutions to the saddle point conditions. They correspond to propagation with 
imaginary momenta, which, as explained above, are neglected in the present level 
of approximation. 

The semiclassical expression for Gq is now substituted in (15) and each of 
the terms in the Born series is evaluated by performing the integrals by the 
stationary phase approximation. A tedious but straight forward algebra yields 
an expression in terms of classical trajectories. To understand the conditions 
which these trajectories should satisfy we should recall that the channel quantum 
numbers m are the quantized classical action variables (measured in units of fi)\ 

j ^ \K\D 

TT 

where ky is the transverse momentum. The conjugate angle variables are 


(63) 


^ Sign(fc^)y7r 
D 

The relevant trajectories emerge from the line a,t x = L with I initial = and 
return to x = L with I final = knk In terms of these trajectories one gets: 









(65) 
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Where the sum is over all the trajectories which satisfy the boundary conditions 
specihed above. is the reduced length : 


— Is k^fi' L (j) sTTi (j) ^TJl 


( 66 ) 


with (ps and cp'^ the initial and hnal angle variables (conjugate to m and m'). Ig is 
the total length of the trajectory s, and Pg is the number of reflections from the 
boundaries. 

Equation (65) plays a central role in the semiclassical discussion. Its right 
hand side looks formally as the standard semiclassical expression for a unitary 
operator. (See e.g. [9]). The unitarity can be checked by performing the matrix 
multiplication by the stationary phase approximation, and therefore the rhs of 
(65)is “semiclassically unitary”. On the Ihs of (65) we hnd S — I which becomes 
unitary only if the identity operator is added to it. The semiclassical approxi¬ 
mation tends to miss this term (compare e.g. the case of elastic scattering from 
a hard disc in the plane. There the semiclassical T matrix is unitary, and the 
missing identity operator (forward scattering amplitude) is missing in the strictly 
semiclassical treatment). In the same spirit we shall omit the identity operator 
from the Ihs of (65). 

Multiplying the semiclassical expression for by we get the semi¬ 

classical expression for S. It coincides with the semiclassical expression for the 
evolution operator corresponding to the mapping of the transverse section at 
X = L on itself by trajectories which start a,t x = L with an initial transverse 
action li and return to x = L after scattering from the boundaries ni x < L 
with transverse action If. To calculate S'” semiclassically, one has to perform 
the intermediate sums by the stationary phase method, and one hnds that the 
stationary point condition corresponds to a reflection of the trajectory at the wall 
X = L. That is, the transverse action variable is unaltered upon reflection, while 
the corresponding angle variable changes its sign. The Trace operation picks up 
periodic orbits, which were reflected from x = L exactly n times. These are 
periodic orbits of the original Sinai billiard, but they are chosen according to the 
multiplicity of their reflections from the transverse line at x = L. We hnally get 
the semiclassical expression for the TrS^ as 


TrS^' 


Ps 


det(/-Mj)|i/2 


r(^lg-\-7ri's ) 


(67) 


Where the sum is over all primitive periodic orbits of the billiard which scatter 
Ps times from x = L and which satisfy n = pgV. Their length is Ig and Mg is the 
monodromy matrix. The index Pg counts the total number of reflections from the 
billiard boundaries (including the n reflections from the line x = L). This ex¬ 
pression completes the semiclassical construction since it can be now substituted 
in (56, 57, 60) to yield the semiclassical expression for the secular equation. 
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Figure 17: Fourier transform of TrS^(k) for (a) = 1 and (b) N = 2 (solid 

line). The shape of the billiard was rendered rectangular (L^ = 2^ Ly = 1, 
R = 0.75) in order to remove degeneracies in the length .spectrum. The dashed 
line .shows the semiclassical prediction obtained from some of the shortest unstable 
PO’s. 


5 Discussion 

The main purpose of this article was to check and demonstrate a few ideas con¬ 
cerning the quantization of billiards using the scattering approach. The system 
which we considered in detail was the Sinai billiard and we addressed it from 
various angles. An important point in the present work was the detailed inves¬ 
tigation of the role of evanescent modes. In this context we would like to recall 
that the Sinai billiard was previously investigated in [5], where it was quantized 
by a similar approach, but for one difference, namely, that the channel was not 
connected to the line x = T but to the section [i?, T] of the line x = 0 (see hgure 
(1)). It is gratifying that the properties of the system are independent of the 
method of attaching the auxiliary continuum (the channel). We would like to 
emphasize that the two scattering problems are very different in character. The 
attachment of the channel to the narrow section (at x = 0) leads to a scattering 
system which is chaotic, since the set of trapped orbits is a non trivial fractal set. 
This is not the case when the channel is connected to the right (at x = T). The 
only trapped orbit in this case is the one which reflects between the circle and the 
straight wall of the billiard along the line x = 0. Hence, the scattering problem 
is not chaotic. Another important difference was explained when the effect of 
evanescent modes was discussed. It was shown that the increase of the radius has 
opposite effects in the two quantization schemes. In the present scheme one had 
to increase the number of evanescent modes as the radius was increased. In the 
method of [5], the effect of the evanescent modes was diminished. 
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The scattering approach developed in this article is related in many ways 
to the method proposed by Bogomolny [8]. As a matter of fact, Bogomolny’s 
approach for billiard problems results in exactly the same semiclassical theory 
as the one obtained here. The main difference is that here we made an explicit 
connection with the scattering problem, thus giving the T operator a physical 
interpretation. The discussion at the end of the previous chapter on the unitarity 
of the semiclassical S matrix, is relevant to Bogomolny’s approach, since also in 
his case, the T matrix is only semiclassically unitarity. The connection with a 
scattering problem sheds light on this phenomenon. 

The scattering approach to quantization, since it starts from an exact expres¬ 
sion can serve as a starting point for deriving improved semiclassical quantization 
procedures. The hrst issue to be addressed is the inclusion of the leading correc¬ 
tions which are coming from evanescent modes. Even though these corrections 
are exponential in h, their effect appear in practical calculations as can be seen 
in hgures (13, 14, 15). Other semiclassical corrections, especially those which 
were recently discussed by Gaspard [20] can be easily incorporated within this 
approach. 

Finally we would like to remark that the same method of quantization which 
was developed here can be easily extended to two dimensional Hamiltonian prob¬ 
lems of the type 

H{q,p) = ^p‘^ + V{^ (68) 

where V{^ is a binding potential so that V^ oo for q —)• oo. This is done 
by dehning an auxiliary scattering problem in the following way. Assume that 
the potential V{^ assumes its minimum value at (f = 0. Construct a scattering 
problem dehned by the potential V = V{^ for qi < 0, and V{^ = V{qi = 
0,q2) for qi > 0. Thus, a channel is obtained, and a scattering matrix can be 
calculated. A corresponding is built by using the same construction as above 
for the right part of the potential V{^. Their product dehnes the S matrix for 
the secular equation det[I — S{E)) = 0. A detailed discussion of this construct 
is presented elsewhere. 
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A Angular momentum decomposition for ex¬ 
ponentials 

Suppose -\- ky = > 0, but k^ and ky not necessarily real. We are looking 

for the radial functions i?/(r) in the angular momentum decomposition of 2d 
exponential functions 


OO 

exp(zA;^a; + ikyy) = ^ Ri{r)e'‘‘^ (69) 

l= — oo 


with r = + j/2, cos ^ = xjr and sin ^ = yjr. The radial function Ri{r) is 

a solution of Bessels differential equation and therefore a linear combination of 
Ji{kr) and Yi{kr) [11]-(9.1). 

Ri{r) = riJi{kr) R siYi{kr) (70) 


Moreover, it is regular in the origin and thus contains no contribution from the 
divergent Neumann function: si = 0. The remaining coefficient r; will be obtained 
from the asymptotic behaviour for r —)• 0. We have in this region [If] 

■Ukr) ^ (/ > 0 ) ( 71 ) 

and 

J.i{kr) = {-ifMkr) (72) 

On the other hand, the Taylor expansion of the exponential reads 


eyr\){ik^^x + ikyy) 


OO 

J^iik^x + ikyijy 

°° {iry 

Y Wi <t^ + ky sin ^)' 

1=0 

tyyYieUe-‘y-,k,{e‘^-e-‘y)' 

1=0 ^ 

°° (jrV / . . X ; 

/=o ^ 


After expanding the sum in powers of and neglecting all but the lowest term 
in r this yields 


exp(zA;^a; + ikyy) Y 


ir 


W2'''K|! 


{kx — i Sign{l)kyy^'^e 


il(f) 


( 74 ) 
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Comparing this to (71) and (72) we obtain 


ri 


i Sign[l)kx + ky 
k 


( 75 ) 


which was the aim of this appendix. If both, k^ and ky^ are real, one can introduce 
(f)k through cos = k^Jk and sin^^ = ky/k and obtains the well known formula 

OO 

eyi\){ikxX + ikyy) = ^ (76) 

l= — oo 


B Symmetries for the generalized S 

In this appendix we show, how unitarity and symmetry of a scattering operator 
can be completed by similiar relations for the nonunitary parts of the generalized 
S. Moreover we will check the implications of these relations for the numerical 
calculation of in the system under consideration. 

In order to get a compact notation we divide the full S into four parts accord¬ 
ing to whether the involved indices correspond to open (o) or closed (c) channels: 


S' 


q q 

^00 

s™ V 


( 77 ) 


Because of the real potential, together with (9) it’s complex conjugate must 
be a solution of our problem as well. We assume an evanescent incoming mode 
n > A and get 




* 

n 


-iknX A -ikmX 

- Y. S*mn fT- 

V m = l V 


OO 

^ /TT— 

n'=A-|-l V ^n' 


( 78 ) 


Any linear combination of and is again a solution of our scattering problem. 
Consider 


A ^-ikmX 

= lY + 

m = l V^m 

A pikmX OO pif^n'X 

+ Y^ / 7 —+ Y^ (Sri'n + S*i^)—i=(/)n' (79) 

m = l V^m ra'=A-|-l VW' 

As for any solution, S will transform the coefficients of the incoming waves into 
those of the outgoing ones. Applying this to the last equation we obtain two 
relations between the different parts of the full S: 

Soc = i SooS*^ (80) 
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and 


2 3(5W-) = 5'.„5': 


( 81 ) 


Doing the same for an open incoming mode one obtains the well known relation 


I 


SooS: 


( 82 ) 


A similiar idea [19] can be used to generalize the unitarity of S. Here we start 
from (33) and use a rectangular area Q enclosed by the channel walls at x = 0, 
y = 0 and y = L together with the line x = Xq. If now and = ^n’ are 

two solutions in the form (9) only the line x = Xq can contribute to the integral. 
Moreover, when integrated over y the double sum of modes will collapse into a 
single sum and we obtain 


0 = + 




^ ^ mn fTT" 


m = l 


+ E 


m = l 




Q ^ ^-iknX 


AknX '' 


cv I ^nn' rj — ^nn' rj — 

UX y ykji 


1 ^mn' 
OX \ 


— tkyj 


ikmX \ 






cv I ^mn 
ox \ 






- s: 


(83) 


There are three different ways to choose n and n': both open, one open and both 
closed. Each will give an additional relation for S after the evaluation of the 
complex conjugate in the last equation. They read as follows: 

- n, n’ open: 


I = SlSl (84) 

Together with (82) this implies symmetry and unitarity for the submatrix of S 
corresponding to open channels. 

- n open, n’ closed: 


Si = z (85) 

Comparing this relation to the previously found (80) we obtain the symmetry for 
the mixed part of S: 

Soc = si ( 86 ) 

- n, n’ closed: 

- Ec = * (87) 
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(88) 


From this we derive using (86) and (81) 

q _ OT 
‘'cc — ‘'cc 

Thus we have established the symmetry of the full S 

S = (89) 

and can use it and the additional relations (80) - (82) for testing the numerical 
accuracy of the calculations. Such a check is particularly important if one deals 
with evanescent modes since the exponential divergence of the involved numbers 
might otherwise lead to large numerical errors. However, the obtained relations 
do not provide a complete check for the generalized S - matrix. To make this 
point clear we consider the restrictions following from (80) - (82) for S in (42). 
For this we divide the matrix (7 in a fashion similiar to (77) into two parts 

C=(^") (90) 

and keep in mind that from the dehnition (38) we have 

c*: = a (91) 

and 

c*: = za. (92) 

Now we hnd from unitarity (82) 

I = I + Co (p-'C-J Co{p-y + p-^ + (P-')*) Cj (93) 

0 = p-^Cj CoiP-y + p-^ + (p-^)* (94) 

0 = Cj Co + P* + P (95) 

»(-?) = -Icj C, (96) 

Similiarly the consistency relation for the mixed part (80) yields 

p-1 Cj = I Co ((P-')* + p-' cj Co [P-y ) {Cjy (97) 

p-^ = -(p-^)* - p-^ C'J^(P”^)* (98) 

and from this follows after multiplication by P and P* once again (96). The 
consistency for the purely closed part (81) results in 

- iCoP-^cj- iiCoyiP-y*icjy = -iCoP-^cyCo{p-y*{cly ( 99 ) 

p-^ + (p-^)* = - p-^ cyco{p~y* (100) 

which is also equivalent to (96). We conclude, that the crucial point in the 
calculation of S is the computation of 71(P), the accuracy of which cannot be 
checked by any of the relations. It contains essentially the imaginary part of the 
structure function and therefore the algorithm described in the appendix of [10] 
for its calculation will determine the accuracy of our eigenvalues. 
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